Fish Comparisons Simple Models

So here we’re going to do some simpler model selection things with the data using GLMs instead of GAMs

Sets up the functions that will be used later

rad2deg <- function(rad) {(rad * 180) / (pi)}
deg2rad <- function(deg) {(deg * pi) / (180)}
round_any <- function(x, accuracy, f=round){f(x/ accuracy) * accuracy}
ang_mean <- function(x){rad2deg(atan(mean(sin(deg2rad(x)))/mean(cos(deg2rad(x)))))}

fold_angle_0_360_to_0_180 <- function(x){abs(abs(x-180)-180)}

Reads in the data and alters it as needed

comp_data <- read.csv("Fish_Comp_Values.csv")
comp_data <- na.omit(comp_data)

comp_data <- comp_data %>% mutate(Flow = ifelse(Flow == "0", "Flow 0", "Flow 2")) %>%
                           mutate(Ablation = ifelse(Ablation == "N", "No Ablation", "Ablated")) %>%
                           mutate(Darkness = ifelse(Darkness == "N", "Light", "Dark")) %>%
                           mutate(Heading_Diff = abs(Heading_Diff)) %>%
                           filter(Distance <= 4) %>%
                           filter(Speed_Diff <= 6) %>% 
                           mutate(quarter_heading_diff = abs(abs(Heading_Diff-90)-90)) %>%
                           mutate(Is_Aligned = ifelse(Heading_Diff < 30, 1, 0)) %>%
                           mutate(Is_Reversed = ifelse(Heading_Diff > 150, 1, 0)) %>%
                           mutate(Angle = fold_angle_0_360_to_0_180(Angle))

sum_comp_data <- comp_data %>% mutate(X_Distance = round_any(X_Distance,0.25), 
                                      Y_Distance = round_any(abs(Y_Distance),0.25)) %>%
                               group_by(Flow,Ablation,Darkness,X_Distance,Y_Distance) %>%
                               summarise(Speed_Diff = mean(Speed_Diff),
                                         Heading_Diff = ang_mean(Heading_Diff),
                                         Sync = mean(Sync),
                                         Quarter_Heading_Diff = mean(quarter_heading_diff),
                                         Is_Aligned = mean(Is_Aligned),
                                         Is_Reversed = mean(Is_Reversed))

`summarise()` has grouped output by 'Flow', 'Ablation', 'Darkness', 'X_Distance'. You can override using the `.groups` argument.

#Basic Stats

speed_anova <- glm(Speed_Diff ~ Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness, data = comp_data)

Anova(speed_anova)

Analysis of Deviance Table (Type II tests)

Response: Speed_Diff
              LR Chisq Df Pr(>Chisq)    
Flow           20.1114  1  7.306e-06 ***
Ablation       24.4652  1  7.566e-07 ***
Darkness       13.7065  1  0.0002137 ***
Flow:Ablation   0.0032  1  0.9545865    
Flow:Darkness   1.5440  1  0.2140271    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

heading_anova <- glm(quarter_heading_diff ~ Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness, data = comp_data)

Anova(heading_anova)

Analysis of Deviance Table (Type II tests)

Response: quarter_heading_diff
              LR Chisq Df Pr(>Chisq)
Flow           0.39037  1     0.5321
Ablation       0.11881  1     0.7303
Darkness       0.44816  1     0.5032
Flow:Ablation  0.44387  1     0.5053
Flow:Darkness  0.32170  1     0.5706

sync_anova <- glm(Sync ~ Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness, data = comp_data)

Anova(sync_anova)

Analysis of Deviance Table (Type II tests)

Response: Sync
              LR Chisq Df Pr(>Chisq)    
Flow            2.0144  1  0.1558116    
Ablation       14.0452  1  0.0001785 ***
Darkness        0.7801  1  0.3771189    
Flow:Ablation   9.3099  1  0.0022792 ** 
Flow:Darkness   7.8318  1  0.0051335 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

dist_anova <- glm(Distance ~ Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness, data = comp_data)

Anova(dist_anova)

Analysis of Deviance Table (Type II tests)

Response: Distance
              LR Chisq Df Pr(>Chisq)    
Flow             0.358  1     0.5496    
Ablation         0.068  1     0.7944    
Darkness       107.389  1     <2e-16 ***
Flow:Ablation    1.970  1     0.1605    
Flow:Darkness    0.032  1     0.8577    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ggplot(comp_data, aes(x = Flow, y = Speed_Diff, fill = paste(Ablation,Darkness,sep="-")))+
  geom_boxplot(outlier.shape = NA) +
  ylim(0,3) +
  theme_light()

ggplot(comp_data, aes(x = Flow, y = quarter_heading_diff, fill = paste(Ablation,Darkness,sep="-")))+
  geom_boxplot(outlier.shape = NA) +
  ylim(0,90) +
  theme_light()

ggplot(comp_data, aes(x = Flow, y = Sync, fill = paste(Ablation,Darkness,sep="-")))+
  geom_boxplot(outlier.shape = NA) +
  ylim(0,1) +
  theme_light()

ggplot(comp_data, aes(x = Flow, y = Distance, fill = paste(Ablation,Darkness,sep="-")))+
  geom_boxplot(outlier.shape = NA) +
  ylim(0,4) +
  theme_light()

ggplot(comp_data, aes(x = Flow, y = Angle, fill = paste(Ablation,Darkness,sep="-")))+
  geom_boxplot(outlier.shape = NA) +
  ylim(0,180) +
  theme_light()

Now we add in distance and angle

speed_anova_da <- glm(Speed_Diff ~ poly(Distance,4)*Angle*(Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness), data = comp_data)

Anova(speed_anova_da)

Analysis of Deviance Table (Type II tests)

Response: Speed_Diff
                                      LR Chisq Df Pr(>Chisq)    
poly(Distance, 4)                      115.177  4  < 2.2e-16 ***
Angle                                    5.946  1  0.0147510 *  
Flow                                    27.012  1  2.022e-07 ***
Ablation                                27.929  1  1.258e-07 ***
Darkness                                 6.729  1  0.0094872 ** 
poly(Distance, 4):Angle                  1.964  4  0.7423022    
Flow:Ablation                            0.006  1  0.9399196    
Flow:Darkness                            0.300  1  0.5840184    
poly(Distance, 4):Flow                   1.872  4  0.7592434    
poly(Distance, 4):Ablation              18.003  4  0.0012324 ** 
poly(Distance, 4):Darkness              20.232  4  0.0004494 ***
Angle:Flow                               5.300  1  0.0213266 *  
Angle:Ablation                           1.855  1  0.1732289    
Angle:Darkness                           0.750  1  0.3863939    
poly(Distance, 4):Flow:Ablation         24.820  4  5.468e-05 ***
poly(Distance, 4):Flow:Darkness          8.276  4  0.0819669 .  
Angle:Flow:Ablation                      0.000  1  0.9967220    
Angle:Flow:Darkness                      0.068  1  0.7937129    
poly(Distance, 4):Angle:Flow            14.614  4  0.0055717 ** 
poly(Distance, 4):Angle:Ablation        12.121  4  0.0164745 *  
poly(Distance, 4):Angle:Darkness         3.761  4  0.4393154    
poly(Distance, 4):Angle:Flow:Ablation   26.352  4  2.687e-05 ***
poly(Distance, 4):Angle:Flow:Darkness    9.819  4  0.0435844 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

heading_anova_da <- glm(quarter_heading_diff ~ Distance*Angle*(Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness), data = comp_data)

Anova(heading_anova_da)

Analysis of Deviance Table (Type II tests)

Response: quarter_heading_diff
                             LR Chisq Df Pr(>Chisq)  
Distance                       2.5860  1    0.10781  
Angle                          0.0237  1    0.87761  
Flow                           0.5212  1    0.47035  
Ablation                       0.1588  1    0.69030  
Darkness                       0.3356  1    0.56238  
Distance:Angle                 2.1602  1    0.14162  
Flow:Ablation                  0.5623  1    0.45332  
Flow:Darkness                  0.1771  1    0.67388  
Distance:Flow                  0.0729  1    0.78715  
Distance:Ablation              1.2093  1    0.27148  
Distance:Darkness              0.1028  1    0.74850  
Angle:Flow                     0.0598  1    0.80688  
Angle:Ablation                 0.9390  1    0.33254  
Angle:Darkness                 0.9958  1    0.31832  
Distance:Flow:Ablation         0.2560  1    0.61286  
Distance:Flow:Darkness         1.4412  1    0.22995  
Angle:Flow:Ablation            0.2370  1    0.62640  
Angle:Flow:Darkness            0.0116  1    0.91409  
Distance:Angle:Flow            0.1515  1    0.69713  
Distance:Angle:Ablation        4.7465  1    0.02936 *
Distance:Angle:Darkness        0.7361  1    0.39092  
Distance:Angle:Flow:Ablation   1.3680  1    0.24216  
Distance:Angle:Flow:Darkness   6.3862  1    0.01150 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

sync_anova_da <- glm(Sync ~ Distance*Angle*(Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness), data = comp_data)

Anova(sync_anova_da)

Analysis of Deviance Table (Type II tests)

Response: Sync
                             LR Chisq Df Pr(>Chisq)    
Distance                       31.817  1  1.694e-08 ***
Angle                           0.666  1  0.4143894    
Flow                            1.399  1  0.2368323    
Ablation                       13.757  1  0.0002081 ***
Darkness                        0.029  1  0.8647383    
Distance:Angle                  3.098  1  0.0783799 .  
Flow:Ablation                   8.824  1  0.0029737 ** 
Flow:Darkness                   7.255  1  0.0070689 ** 
Distance:Flow                   1.241  1  0.2652318    
Distance:Ablation               0.036  1  0.8487728    
Distance:Darkness               1.186  1  0.2761968    
Angle:Flow                      2.525  1  0.1120840    
Angle:Ablation                  5.473  1  0.0193101 *  
Angle:Darkness                  1.156  1  0.2822980    
Distance:Flow:Ablation          7.958  1  0.0047885 ** 
Distance:Flow:Darkness          2.776  1  0.0956635 .  
Angle:Flow:Ablation             1.892  1  0.1689323    
Angle:Flow:Darkness             0.068  1  0.7941314    
Distance:Angle:Flow             0.201  1  0.6536996    
Distance:Angle:Ablation        11.861  1  0.0005733 ***
Distance:Angle:Darkness         0.150  1  0.6983841    
Distance:Angle:Flow:Ablation    0.871  1  0.3506822    
Distance:Angle:Flow:Darkness    7.252  1  0.0070829 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ggplot(comp_data, aes(x = Distance, y = Speed_Diff, color = paste(Ablation,Darkness,sep="-"), fill = paste(Ablation,Darkness,sep="-")))+
  geom_smooth(method = "glm", formula = y ~ poly(x,4))+
  facet_wrap(~ Flow) +
  theme_light()

ggplot(comp_data, aes(x = Distance, y = quarter_heading_diff, color = paste(Ablation,Darkness,sep="-"), fill = paste(Ablation,Darkness,sep="-")))+
  geom_smooth(method = "glm", formula = y ~ poly(x,4))+
  facet_wrap(~ Flow) +
  theme_light()

ggplot(comp_data, aes(x = Distance, y = Sync, color = paste(Ablation,Darkness,sep="-"), fill = paste(Ablation,Darkness,sep="-")))+
  geom_smooth(method = "glm", formula = y ~ poly(x,4))+
  facet_wrap(~ Flow) +
  theme_light()

ggplot(comp_data, aes(x = Angle, y = Speed_Diff, color = paste(Ablation,Darkness,sep="-"), fill = paste(Ablation,Darkness,sep="-")))+
  geom_smooth(method = "glm", formula = y ~ poly(x,4))+
  facet_wrap(~ Flow) +
  theme_light()

ggplot(comp_data, aes(x = Angle, y = quarter_heading_diff, color = paste(Ablation,Darkness,sep="-"), fill = paste(Ablation,Darkness,sep="-")))+
  geom_smooth(method = "glm", formula = y ~ poly(x,4))+
  facet_wrap(~ Flow) +
  theme_light()

ggplot(comp_data, aes(x = Angle, y = Sync, color = paste(Ablation,Darkness,sep="-"), fill = paste(Ablation,Darkness,sep="-")))+
  geom_smooth(method = "glm", formula = y ~ poly(x,4))+
  facet_wrap(~ Flow) +
  theme_light()

This creates a dataframe that is used for making model predictions and graphing them

d <- seq(from = 0, to = 4, by = 0.1)
a <- seq(from = 0, to = 180, by = 10)

flows <- c("Flow 0", "Flow 2")
ablation <- c("No Ablation", "Ablated")
dark <- c("Light","Dark")

predict_df_da <- expand.grid(Distance = d, Angle = a, Flow = flows, Ablation = ablation, Darkness = dark)
predict_df_da <- predict_df_da %>% mutate(X_Distance = Distance*(cos(deg2rad(Angle))), Y_Distance = Distance*(sin(deg2rad(Angle))))

predict_df_da <- predict_df_da %>% filter(!(Ablation == "Ablated" & Darkness == 'Dark'))
predict_df_da <- na.omit(predict_df_da)

Speed Models

So I want to make a model that tries to include both distance and angle. What I really need is to figure out what power works best for angle and distance, and then include that. The issue is stepwise selection can just remove one of them entirely. So I’m going to use cross validation to select the power for both of them, and then use that. I will also calculate the AICs to help in this decision making process

set.seed(7)

comp_data_speed_model <- comp_data %>% select(c(Flow,Darkness,Ablation,Angle,Distance,Speed_Diff))

max_poly <- 10

speed_dist_cv_error_10 <- rep(0,max_poly)
speed_dist_AIC_10 <- rep(0,max_poly)


for (i in 1:max_poly){
  speed_dist_fit <- glm(Speed_Diff ~ poly(Distance,i)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
                        data = comp_data_speed_model)
  speed_dist_cv_error_10[i] <- cv.glm(comp_data_speed_model, speed_dist_fit, K = 10)$delta[1]
  speed_dist_AIC_10[i] <- AIC(speed_dist_fit)
}

speed_angle_cv_error_10 <- rep(0,max_poly)
speed_angle_AIC_10 <- rep(0,max_poly)


for (i in 1:max_poly){
  speed_angle_fit <- glm(Speed_Diff ~ poly(Angle,i)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
                         data = comp_data_speed_model)
  speed_angle_cv_error_10[i] <- cv.glm(comp_data_speed_model, speed_angle_fit, K = 10)$delta[1]
  speed_angle_AIC_10[i] <- AIC(speed_angle_fit)
}

speed_poly_plot_df <- data.frame(c(seq(max_poly),seq(max_poly)),
                                 c(speed_dist_cv_error_10,speed_angle_cv_error_10),
                                 c(speed_dist_AIC_10,speed_angle_AIC_10),
                                 c(rep("Distance",max_poly),rep("Angle",max_poly)))

colnames(speed_poly_plot_df) <- c("Degree","Error","AIC","Predictor")

speed_poly_plot_df <- speed_poly_plot_df %>% group_by(Predictor) %>%
                                             mutate(minError = min(Error),minAIC = min(AIC)) %>%
                                             ungroup() %>%
                                             mutate(isMinEror = ifelse(Error == minError,3,1),isMinAIC = ifelse(AIC == minAIC,3,1))

ggplot(speed_poly_plot_df, aes(x = Degree, y = Error, color = Predictor))+
  geom_point(size = speed_poly_plot_df$isMinEror)+
  geom_line()+
  theme_light()+
  facet_wrap(~ Predictor, scales = "free") +
  scale_size(guide = "none")

ggplot(speed_poly_plot_df, aes(x = Degree, y = AIC, color = Predictor))+
  geom_point(size = speed_poly_plot_df$isMinAIC)+
  geom_line()+
  theme_light()+
  facet_wrap(~ Predictor, scales = "free") +
  scale_size(guide = "none")

Now let’s use those polynomials to create the actual model

While they not be the lowest, the big improvements happen at 4 for Distance, and 2 for Angle

speed_diff_glm <- glm(Speed_Diff ~ poly(Distance,4)*poly(Angle,2)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
                        data = comp_data_speed_model)

summary(speed_diff_glm)

Call:
glm(formula = Speed_Diff ~ poly(Distance, 4) * poly(Angle, 2) * 
    (Ablation + Flow + Darkness + Ablation:Flow + Darkness:Flow), 
    data = comp_data_speed_model)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.2116  -0.3067  -0.0859   0.1910   3.4546  

Coefficients:
                                                                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)                                                        8.374e-01  3.163e-02  26.476  < 2e-16 ***
poly(Distance, 4)1                                                 5.475e+00  3.118e+00   1.756 0.079188 .  
poly(Distance, 4)2                                                 6.853e+00  3.187e+00   2.150 0.031583 *  
poly(Distance, 4)3                                                 1.147e+00  2.850e+00   0.402 0.687379    
poly(Distance, 4)4                                                 7.173e+00  2.413e+00   2.973 0.002962 ** 
poly(Angle, 2)1                                                    4.186e+00  2.370e+00   1.766 0.077374 .  
poly(Angle, 2)2                                                   -2.217e+00  2.752e+00  -0.806 0.420391    
AblationNo Ablation                                                7.305e-02  2.163e-02   3.377 0.000736 ***
FlowFlow 2                                                         5.948e-02  4.749e-02   1.253 0.210403    
DarknessLight                                                     -4.160e-02  2.755e-02  -1.510 0.131132    
poly(Distance, 4)1:poly(Angle, 2)1                                 4.232e+02  2.309e+02   1.832 0.066938 .  
poly(Distance, 4)2:poly(Angle, 2)1                                 3.148e+02  2.369e+02   1.329 0.183886    
poly(Distance, 4)3:poly(Angle, 2)1                                 3.883e+02  2.266e+02   1.713 0.086702 .  
poly(Distance, 4)4:poly(Angle, 2)1                                -9.528e+01  2.035e+02  -0.468 0.639653    
poly(Distance, 4)1:poly(Angle, 2)2                                -1.237e+02  3.028e+02  -0.409 0.682821    
poly(Distance, 4)2:poly(Angle, 2)2                                 2.186e+02  3.088e+02   0.708 0.479134    
poly(Distance, 4)3:poly(Angle, 2)2                                 3.175e+02  2.579e+02   1.231 0.218397    
poly(Distance, 4)4:poly(Angle, 2)2                                 8.375e+01  2.127e+02   0.394 0.693829    
AblationNo Ablation:FlowFlow 2                                     1.443e-02  3.753e-02   0.384 0.700689    
FlowFlow 2:DarknessLight                                           4.637e-03  4.152e-02   0.112 0.911073    
poly(Distance, 4)1:AblationNo Ablation                            -5.285e+00  2.275e+00  -2.323 0.020220 *  
poly(Distance, 4)2:AblationNo Ablation                            -2.449e+00  2.303e+00  -1.063 0.287643    
poly(Distance, 4)3:AblationNo Ablation                            -1.821e+00  2.002e+00  -0.910 0.362897    
poly(Distance, 4)4:AblationNo Ablation                            -3.468e+00  1.604e+00  -2.162 0.030663 *  
poly(Distance, 4)1:FlowFlow 2                                     -7.224e+00  4.631e+00  -1.560 0.118819    
poly(Distance, 4)2:FlowFlow 2                                     -1.121e+01  4.809e+00  -2.330 0.019838 *  
poly(Distance, 4)3:FlowFlow 2                                     -4.884e-01  4.382e+00  -0.111 0.911259    
poly(Distance, 4)4:FlowFlow 2                                     -4.360e-01  3.740e+00  -0.117 0.907201    
poly(Distance, 4)1:DarknessLight                                   2.302e+00  2.659e+00   0.866 0.386686    
poly(Distance, 4)2:DarknessLight                                  -5.323e+00  2.769e+00  -1.922 0.054593 .  
poly(Distance, 4)3:DarknessLight                                  -1.956e+00  2.514e+00  -0.778 0.436673    
poly(Distance, 4)4:DarknessLight                                  -4.435e+00  2.144e+00  -2.069 0.038619 *  
poly(Angle, 2)1:AblationNo Ablation                               -3.353e+00  1.624e+00  -2.065 0.039006 *  
poly(Angle, 2)2:AblationNo Ablation                                4.392e+00  1.976e+00   2.223 0.026278 *  
poly(Angle, 2)1:FlowFlow 2                                        -1.248e+00  3.654e+00  -0.342 0.732653    
poly(Angle, 2)2:FlowFlow 2                                         8.118e+00  4.250e+00   1.910 0.056166 .  
poly(Angle, 2)1:DarknessLight                                     -1.067e+00  2.076e+00  -0.514 0.607152    
poly(Angle, 2)2:DarknessLight                                     -6.070e-01  2.352e+00  -0.258 0.796389    
poly(Distance, 4)1:AblationNo Ablation:FlowFlow 2                  1.004e+01  3.834e+00   2.617 0.008880 ** 
poly(Distance, 4)2:AblationNo Ablation:FlowFlow 2                  4.709e+00  3.926e+00   1.199 0.230398    
poly(Distance, 4)3:AblationNo Ablation:FlowFlow 2                  1.692e-02  3.463e+00   0.005 0.996103    
poly(Distance, 4)4:AblationNo Ablation:FlowFlow 2                 -2.326e-01  2.877e+00  -0.081 0.935552    
poly(Distance, 4)1:FlowFlow 2:DarknessLight                        4.597e+00  3.940e+00   1.167 0.243451    
poly(Distance, 4)2:FlowFlow 2:DarknessLight                        1.025e+01  4.104e+00   2.498 0.012525 *  
poly(Distance, 4)3:FlowFlow 2:DarknessLight                       -5.024e-01  3.813e+00  -0.132 0.895185    
poly(Distance, 4)4:FlowFlow 2:DarknessLight                        5.413e-01  3.249e+00   0.167 0.867682    
poly(Angle, 2)1:AblationNo Ablation:FlowFlow 2                     1.924e+00  2.904e+00   0.663 0.507641    
poly(Angle, 2)2:AblationNo Ablation:FlowFlow 2                    -6.275e+00  3.508e+00  -1.789 0.073740 .  
poly(Angle, 2)1:FlowFlow 2:DarknessLight                           8.569e-01  3.212e+00   0.267 0.789669    
poly(Angle, 2)2:FlowFlow 2:DarknessLight                           1.495e+00  3.659e+00   0.409 0.682838    
poly(Distance, 4)1:poly(Angle, 2)1:AblationNo Ablation            -3.168e+02  1.646e+02  -1.925 0.054294 .  
poly(Distance, 4)2:poly(Angle, 2)1:AblationNo Ablation            -4.343e+02  1.636e+02  -2.655 0.007962 ** 
poly(Distance, 4)3:poly(Angle, 2)1:AblationNo Ablation            -9.630e+01  1.515e+02  -0.636 0.525092    
poly(Distance, 4)4:poly(Angle, 2)1:AblationNo Ablation             5.509e+01  1.300e+02   0.424 0.671750    
poly(Distance, 4)1:poly(Angle, 2)2:AblationNo Ablation             1.351e+02  2.312e+02   0.584 0.559017    
poly(Distance, 4)2:poly(Angle, 2)2:AblationNo Ablation            -5.181e+01  2.345e+02  -0.221 0.825118    
poly(Distance, 4)3:poly(Angle, 2)2:AblationNo Ablation            -1.465e+02  1.912e+02  -0.766 0.443742    
poly(Distance, 4)4:poly(Angle, 2)2:AblationNo Ablation             1.054e+02  1.488e+02   0.708 0.478770    
poly(Distance, 4)1:poly(Angle, 2)1:FlowFlow 2                     -7.079e+02  3.478e+02  -2.035 0.041884 *  
poly(Distance, 4)2:poly(Angle, 2)1:FlowFlow 2                     -9.543e+01  3.550e+02  -0.269 0.788075    
poly(Distance, 4)3:poly(Angle, 2)1:FlowFlow 2                     -5.019e+02  3.380e+02  -1.485 0.137620    
poly(Distance, 4)4:poly(Angle, 2)1:FlowFlow 2                      1.022e+03  3.057e+02   3.344 0.000832 ***
poly(Distance, 4)1:poly(Angle, 2)2:FlowFlow 2                      3.442e+02  4.573e+02   0.753 0.451688    
poly(Distance, 4)2:poly(Angle, 2)2:FlowFlow 2                     -4.818e+02  4.694e+02  -1.026 0.304771    
poly(Distance, 4)3:poly(Angle, 2)2:FlowFlow 2                     -3.571e+02  4.016e+02  -0.889 0.373974    
poly(Distance, 4)4:poly(Angle, 2)2:FlowFlow 2                     -3.991e+02  3.306e+02  -1.207 0.227334    
poly(Distance, 4)1:poly(Angle, 2)1:DarknessLight                  -7.976e+01  1.999e+02  -0.399 0.689922    
poly(Distance, 4)2:poly(Angle, 2)1:DarknessLight                   1.168e+02  2.100e+02   0.556 0.578190    
poly(Distance, 4)3:poly(Angle, 2)1:DarknessLight                  -2.198e+02  2.025e+02  -1.086 0.277713    
poly(Distance, 4)4:poly(Angle, 2)1:DarknessLight                   1.097e+02  1.827e+02   0.601 0.548090    
poly(Distance, 4)1:poly(Angle, 2)2:DarknessLight                   1.008e+02  2.534e+02   0.398 0.690692    
poly(Distance, 4)2:poly(Angle, 2)2:DarknessLight                  -3.290e+02  2.654e+02  -1.239 0.215268    
poly(Distance, 4)3:poly(Angle, 2)2:DarknessLight                  -1.985e+02  2.247e+02  -0.883 0.377124    
poly(Distance, 4)4:poly(Angle, 2)2:DarknessLight                  -2.237e+02  1.881e+02  -1.189 0.234394    
poly(Distance, 4)1:poly(Angle, 2)1:AblationNo Ablation:FlowFlow 2  2.673e+02  2.841e+02   0.941 0.346800    
poly(Distance, 4)2:poly(Angle, 2)1:AblationNo Ablation:FlowFlow 2  6.864e+01  2.836e+02   0.242 0.808736    
poly(Distance, 4)3:poly(Angle, 2)1:AblationNo Ablation:FlowFlow 2 -6.060e+01  2.623e+02  -0.231 0.817321    
poly(Distance, 4)4:poly(Angle, 2)1:AblationNo Ablation:FlowFlow 2 -9.584e+02  2.302e+02  -4.163 3.18e-05 ***
poly(Distance, 4)1:poly(Angle, 2)2:AblationNo Ablation:FlowFlow 2 -3.361e+02  3.930e+02  -0.855 0.392475    
poly(Distance, 4)2:poly(Angle, 2)2:AblationNo Ablation:FlowFlow 2  1.646e+02  3.993e+02   0.412 0.680267    
poly(Distance, 4)3:poly(Angle, 2)2:AblationNo Ablation:FlowFlow 2  7.080e+01  3.353e+02   0.211 0.832742    
poly(Distance, 4)4:poly(Angle, 2)2:AblationNo Ablation:FlowFlow 2  1.890e+02  2.615e+02   0.723 0.469682    
poly(Distance, 4)1:poly(Angle, 2)1:FlowFlow 2:DarknessLight        3.462e+02  2.970e+02   1.166 0.243717    
poly(Distance, 4)2:poly(Angle, 2)1:FlowFlow 2:DarknessLight       -1.392e+02  3.029e+02  -0.460 0.645818    
poly(Distance, 4)3:poly(Angle, 2)1:FlowFlow 2:DarknessLight        4.288e+02  2.946e+02   1.456 0.145562    
poly(Distance, 4)4:poly(Angle, 2)1:FlowFlow 2:DarknessLight       -6.988e+02  2.675e+02  -2.612 0.009022 ** 
poly(Distance, 4)1:poly(Angle, 2)2:FlowFlow 2:DarknessLight        1.313e+02  3.807e+02   0.345 0.730104    
poly(Distance, 4)2:poly(Angle, 2)2:FlowFlow 2:DarknessLight        4.487e+02  3.929e+02   1.142 0.253424    
poly(Distance, 4)3:poly(Angle, 2)2:FlowFlow 2:DarknessLight        2.905e+02  3.426e+02   0.848 0.396498    
poly(Distance, 4)4:poly(Angle, 2)2:FlowFlow 2:DarknessLight        5.870e+02  2.864e+02   2.050 0.040424 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.2194725)

    Null deviance: 1416.0  on 6050  degrees of freedom
Residual deviance: 1308.3  on 5961  degrees of freedom
AIC: 8086.8

Number of Fisher Scoring iterations: 2

Now let’s make some predictions

speed_pred <- predict_df_da %>% mutate(Speed_Diff = predict(speed_diff_glm,predict_df_da))

comp_data <- comp_data %>% mutate(Round_Dist = as.factor(round_any(Distance,1)), Round_Angle = as.factor(round_any(Angle,30)))

ggplot()+
  geom_boxplot(data = comp_data, aes(x = Round_Dist, y = Speed_Diff))+
  facet_wrap(~ Flow + Ablation + Darkness) +
  theme_light()

round_dist_aov <- aov(Speed_Diff ~ Round_Dist*(Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness), data = comp_data)
Anova(round_dist_aov)

Anova Table (Type II tests)

Response: Speed_Diff
                          Sum Sq   Df F value    Pr(>F)    
Round_Dist                 21.92    4 24.3669 < 2.2e-16 ***
Flow                        5.30    1 23.5595 1.242e-06 ***
Ablation                    5.98    1 26.5765 2.614e-07 ***
Darkness                    1.68    1  7.4726  0.006283 ** 
Flow:Ablation               0.00    1  0.0057  0.939866    
Flow:Darkness               0.10    1  0.4649  0.495347    
Round_Dist:Flow             0.96    4  1.0642  0.372465    
Round_Dist:Ablation         2.19    4  2.4329  0.045321 *  
Round_Dist:Darkness         2.99    4  3.3206  0.010029 *  
Round_Dist:Flow:Ablation    5.31    4  5.9036 9.730e-05 ***
Round_Dist:Flow:Darkness    1.94    4  2.1605  0.070840 .  
Residuals                1354.08 6021                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ggplot()+
  geom_boxplot(data = comp_data, aes(x = Round_Angle, y = Speed_Diff))+
  facet_wrap(~ Flow + Ablation + Darkness) +
  theme_light()

round_angle_aov <- aov(Speed_Diff ~ Round_Angle*(Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness), data = comp_data)
Anova(round_angle_aov)

Anova Table (Type II tests)

Response: Speed_Diff
                           Sum Sq   Df F value    Pr(>F)    
Round_Angle                  3.25    6  2.3765 0.0269962 *  
Flow                         4.90    1 21.5029 3.608e-06 ***
Ablation                     6.21    1 27.2924 1.808e-07 ***
Darkness                     2.99    1 13.1415 0.0002912 ***
Flow:Ablation                0.05    1  0.2310 0.6307756    
Flow:Darkness                0.36    1  1.5723 0.2099197    
Round_Angle:Flow             7.43    6  5.4377 1.280e-05 ***
Round_Angle:Ablation         2.02    6  1.4765 0.1818843    
Round_Angle:Darkness         1.43    6  1.0461 0.3930796    
Round_Angle:Flow:Ablation    5.99    6  4.3881 0.0001975 ***
Round_Angle:Flow:Darkness    0.98    6  0.7154 0.6372082    
Residuals                 1368.14 6009                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# ggplot()+
#   geom_point(data = comp_data %>% filter(Speed_Diff <= 2), aes(x = Distance, y = Speed_Diff), alpha = 0.1)+
#   geom_smooth(data = speed_pred, aes(x = Distance, y = Speed_Diff))+
#   facet_wrap(~ Flow + Ablation + Darkness) +
#   theme_light()

# ggplot()+
#   geom_density_2d_filled(data = comp_data %>% filter(Speed_Diff <= 2), aes(x = Distance, y = Speed_Diff), contour_var = "ndensity")+
#   geom_smooth(data = speed_pred, aes(x = Distance, y = Speed_Diff, color = "red"))+
#   facet_wrap(~ Flow + Ablation + Darkness) +
#   theme_light()

# ggplot()+
#   geom_point(data = comp_data %>% filter(Speed_Diff <= 2), aes(x = Angle, y = Speed_Diff), alpha = 0.1)+
#   geom_smooth(data = speed_pred, aes(x = Angle, y = Speed_Diff))+
#   facet_wrap(~ Flow + Ablation + Darkness) +
#   theme_light()

# ggplot()+
#   geom_density_2d_filled(data = comp_data %>% filter(Speed_Diff <= 2), aes(x = Angle, y = Speed_Diff), contour_var = "ndensity")+
#   geom_smooth(data = speed_pred, aes(x = Angle, y = Speed_Diff, color = "red"))+
#   facet_wrap(~ Flow + Ablation + Darkness) +
#   theme_light()

ggplot(data = comp_data, aes(x = Angle, y = Distance, z = Speed_Diff))+
  stat_summary_2d() +
  facet_wrap(~ Flow + Ablation + Darkness) +
  scale_fill_viridis(direction = -1) +
  theme_light()

You know let’s just try stepwise as well. There’s jsut way too many thing in that possible model

speed_m_all <- glm(Speed_Diff ~ Distance*Angle*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow) +
                                I(Distance^2)*I(Angle^2)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow)+
                                I(Distance^3)*I(Angle^3)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow)+
                                I(Distance^4)*I(Angle^4)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
                        data = comp_data_speed_model)

speed_m_none <- glm(Speed_Diff ~ 1, data = comp_data_speed_model)

speed_m_both <- step(speed_m_none, direction = "both", scope = formula(speed_m_all), trace = F)

summary(speed_m_both)

Call:
glm(formula = Speed_Diff ~ Flow + Distance + I(Angle^4) + Ablation + 
    Angle + I(Angle^3) + I(Angle^2) + Flow:I(Angle^4) + Distance:Ablation + 
    Flow:Ablation + Ablation:Angle + Flow:Angle + Flow:Distance + 
    Flow:Distance:Ablation, data = comp_data_speed_model)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8986  -0.3680  -0.1121   0.2432   3.7254  

Coefficients:
                                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                              1.535e-01  6.916e-02   2.220 0.026451 *  
FlowFlow 2                               2.876e-01  6.243e-02   4.608 4.16e-06 ***
Distance                                 1.190e-01  1.528e-02   7.791 7.77e-15 ***
I(Angle^4)                              -4.649e-09  1.644e-09  -2.828 0.004706 ** 
AblationNo Ablation                      3.053e-01  5.171e-02   5.905 3.73e-09 ***
Angle                                    8.898e-03  3.795e-03   2.345 0.019083 *  
I(Angle^3)                               1.682e-06  6.132e-07   2.743 0.006113 ** 
I(Angle^2)                              -1.965e-04  7.740e-05  -2.539 0.011145 *  
FlowFlow 2:I(Angle^4)                    4.946e-10  1.137e-10   4.350 1.38e-05 ***
Distance:AblationNo Ablation            -9.326e-02  1.966e-02  -4.743 2.15e-06 ***
FlowFlow 2:AblationNo Ablation          -1.349e-01  6.483e-02  -2.081 0.037497 *  
AblationNo Ablation:Angle               -1.082e-03  3.179e-04  -3.404 0.000668 ***
FlowFlow 2:Angle                        -2.065e-03  6.409e-04  -3.222 0.001280 ** 
FlowFlow 2:Distance                     -4.181e-02  2.399e-02  -1.743 0.081447 .  
FlowFlow 2:Distance:AblationNo Ablation  1.065e-01  3.036e-02   3.508 0.000456 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.2679219)

    Null deviance: 1716.9  on 6103  degrees of freedom
Residual deviance: 1631.4  on 6089  degrees of freedom
AIC: 9300

Number of Fisher Scoring iterations: 2

Sync Models

Now let’s try the same thing for sync values (Rayleigh’s R)

set.seed(7)

comp_data_sync_model <- comp_data %>% select(c(Flow,Darkness,Ablation,Angle,Distance,Sync))

max_poly <- 10

sync_dist_cv_error_10 <- rep(0,max_poly)
sybc_dist_AIC_10 <- rep(0,max_poly)


for (i in 1:max_poly){
  sync_dist_fit <- glm(Sync ~ poly(Distance,i)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
                        data = comp_data_sync_model)
  sync_dist_cv_error_10[i] <- cv.glm(comp_data_sync_model, sync_dist_fit, K = 10)$delta[1]
  sybc_dist_AIC_10[i] <- AIC(sync_dist_fit)
}

sync_angle_cv_error_10 <- rep(0,max_poly)
sync_angle_AIC_10 <- rep(0,max_poly)


for (i in 1:max_poly){
  sync_angle_fit <- glm(Sync ~ poly(Angle,i)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
                         data = comp_data_sync_model)
  sync_angle_cv_error_10[i] <- cv.glm(comp_data_sync_model, sync_angle_fit, K = 10)$delta[1]
  sync_angle_AIC_10[i] <- AIC(sync_angle_fit)
}

sync_poly_plot_df <- data.frame(c(seq(max_poly),seq(max_poly)),
                                 c(speed_dist_cv_error_10,speed_angle_cv_error_10),
                                 c(speed_dist_AIC_10,speed_angle_AIC_10),
                                 c(rep("Distance",max_poly),rep("Angle",max_poly)))

colnames(sync_poly_plot_df) <- c("Degree","Error","AIC","Predictor")

sync_poly_plot_df <- sync_poly_plot_df %>% group_by(Predictor) %>%
                                             mutate(minError = min(Error),minAIC = min(AIC)) %>%
                                             ungroup() %>%
                                             mutate(isMinEror = ifelse(Error == minError,3,1),isMinAIC = ifelse(AIC == minAIC,3,1))

ggplot(sync_poly_plot_df, aes(x = Degree, y = Error, color = Predictor))+
  geom_point(size = sync_poly_plot_df$isMinEror)+
  geom_line()+
  theme_light()+
  facet_wrap(~ Predictor, scales = "free") +
  scale_size(guide = "none")

ggplot(sync_poly_plot_df, aes(x = Degree, y = AIC, color = Predictor))+
  geom_point(size = sync_poly_plot_df$isMinAIC)+
  geom_line()+
  theme_light()+
  facet_wrap(~ Predictor, scales = "free") +
  scale_size(guide = "none")

Now let’s use those polynomials to create the actual model

While they not be the lowest, the big improvements happen at 2 for Distance, and 2 for Angle

sync_diff_glm <- glm(Sync ~ I(Distance^2)*I(Angle^2)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
                        data = comp_data_sync_model)

summary(sync_diff_glm)

Call:
glm(formula = Sync ~ I(Distance^2) * I(Angle^2) * (Ablation + 
    Flow + Darkness + Ablation:Flow + Darkness:Flow), data = comp_data_sync_model)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.55625  -0.19623  -0.01846   0.19226   0.51838  

Coefficients:
                                                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                                              4.374e-01  3.651e-02  11.983  < 2e-16 ***
I(Distance^2)                                            8.645e-03  6.445e-03   1.341 0.179875    
I(Angle^2)                                               5.721e-06  2.576e-06   2.221 0.026404 *  
AblationNo Ablation                                      8.300e-02  2.232e-02   3.718 0.000202 ***
FlowFlow 2                                               2.099e-01  5.287e-02   3.971 7.24e-05 ***
DarknessLight                                            5.984e-02  3.305e-02   1.811 0.070222 .  
I(Distance^2):I(Angle^2)                                -7.040e-07  5.514e-07  -1.277 0.201717    
AblationNo Ablation:FlowFlow 2                          -1.648e-01  3.754e-02  -4.389 1.16e-05 ***
FlowFlow 2:DarknessLight                                -9.831e-02  4.694e-02  -2.094 0.036263 *  
I(Distance^2):AblationNo Ablation                       -1.454e-03  4.438e-03  -0.328 0.743149    
I(Distance^2):FlowFlow 2                                -3.488e-02  9.445e-03  -3.692 0.000224 ***
I(Distance^2):DarknessLight                             -1.244e-02  5.622e-03  -2.214 0.026893 *  
I(Angle^2):AblationNo Ablation                          -3.000e-06  1.686e-06  -1.780 0.075200 .  
I(Angle^2):FlowFlow 2                                   -1.480e-05  3.804e-06  -3.891 0.000101 ***
I(Angle^2):DarknessLight                                -3.767e-06  2.288e-06  -1.647 0.099688 .  
I(Distance^2):AblationNo Ablation:FlowFlow 2             2.164e-02  7.423e-03   2.915 0.003567 ** 
I(Distance^2):FlowFlow 2:DarknessLight                   2.284e-02  8.053e-03   2.836 0.004585 ** 
I(Angle^2):AblationNo Ablation:FlowFlow 2                8.795e-06  2.705e-06   3.251 0.001157 ** 
I(Angle^2):FlowFlow 2:DarknessLight                      7.684e-06  3.377e-06   2.275 0.022920 *  
I(Distance^2):I(Angle^2):AblationNo Ablation            -4.831e-08  4.062e-07  -0.119 0.905336    
I(Distance^2):I(Angle^2):FlowFlow 2                      3.123e-06  8.161e-07   3.826 0.000131 ***
I(Distance^2):I(Angle^2):DarknessLight                   7.798e-07  4.726e-07   1.650 0.098985 .  
I(Distance^2):I(Angle^2):AblationNo Ablation:FlowFlow 2 -1.543e-06  6.555e-07  -2.354 0.018600 *  
I(Distance^2):I(Angle^2):FlowFlow 2:DarknessLight       -2.089e-06  7.046e-07  -2.965 0.003042 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.06028771)

    Null deviance: 371.77  on 6103  degrees of freedom
Residual deviance: 366.55  on 6080  degrees of freedom
AIC: 204.5

Number of Fisher Scoring iterations: 2

Now let’s make some predictions

sync_pred <- predict_df_da %>% mutate(Sync = predict(sync_diff_glm,predict_df_da))

ggplot(comp_data, aes(x = Distance, y = Sync))+
  geom_point(alpha = 0.1)+
  geom_smooth(method = lm, formula = y ~ poly(x, 2)) +
  facet_wrap(~ Flow + Ablation + Darkness) +
  theme_light()

ggplot(comp_data, aes(x = Angle, y = Sync))+
  geom_point(alpha = 0.1)+
  geom_smooth(method = glm, formula = y ~ poly(x, 2))+
  facet_wrap(~ Flow + Ablation + Darkness) +
  theme_light()

ggplot(sync_pred, aes(x = Distance, y = Sync))+
  geom_smooth()+
  scale_fill_viridis() +
  facet_wrap(~ Flow + Ablation + Darkness) +
  theme_light()

ggplot(sync_pred, aes(x = Angle, y = Sync))+
  geom_smooth()+
  scale_fill_viridis() +
  facet_wrap(~ Flow + Ablation + Darkness) +
  theme_light()